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We present numerical results demonstrating the possibility of thermalization of single-particle 
observables in a one-dimensional system, which is integrable in both the quantum and classical 
(mean-field) descriptions (a quasicondensate of ultracold, weakly interacting bosonic atoms are 
studied as a definite example). We find that certain initial conditions admit the relaxation of 
single-particle observables to the equilibrium state reasonably close to that corresponding to the 
Bose-Einstein thermal distribution of Bogoliubov quasiparticles. 



PACS numbers: 03.75.Kk,03.75.Gg,02.30.Ik,67.85.-d 



I. INTRODUCTION 



A one-dimensional (ID) system of identical bosons 
with contact interactions is known to be integrable since 
Licb and Liniger have solved analytically the correspond- 
ing quantum problem by means of Bethe ansatz 
In the weakly interacting limit, this system can be de- 
scribed in the mean-field approximation by the Gross- 
Pitaevskii equation (GPE), also known as the nonlinear 
Schrodinger equation (NLSE). Zakharov and Shabat Q 
have demonstrated that the NLSE with defocusing non- 
linearity (which corresponds to the repulsive interactions 
between particles) is integrable by the inverse scattering 
transform (see Q for a general review of the inverse scat- 
tering transform method) . Since the number of integrals 
of motion in an integrable system equals to the number 
of degrees of freedom (infinite in the continuous mean- 
field description @ or equal to the number of particles in 
the quantum Lieb-Linigcr model [![), one might expect 
that the finally attained equilibrium state must still bear 
signatures of the initial conditions. 

One-dimensional bosonic systems have been experi- 
mentally implemented with ultra-cold atoms on atom 
chips 0, H|, with the radial trapping frequency being 
~ 10 3 times higher than the longitudinal one. The ultra- 
cold degenerate atomic system (quasicondensate, i.e. a 
system describable by a macroscopic wave function with 
a fluctuating local phase) was in the ID regime since 
both the temperature and the mean interaction energy 
per atom were well below the energy interval between 
the ground and the first excited states of the radial mo- 
tion. The fact that the static and dynamic correlation 
properties of these systems were in a very good agree- 
ment with the Bose-Einstein equilibrium distribution of 
quasiparticles seemed to be in contradiction with the sys- 
tem intcgrability and called for explanation. To explain 
the observed relaxation of single- and two-particle distri- 
bution functions for the elementary excitations (Bogoli- 
ubov quasiparticles) to the Bose-Einstein equilibrium, a 
mechanism of intcgrability breakdown via three-body ef- 
fective collisions involving virtual excitations of the radial 
degrees of freedom has been proposed Q. 



In the present paper we numerically show the existence 
of a certain case of nonequilibrium initial conditions of 
the GPE, which provide a very fast relaxation of the sim- 
plest (single-particle) observables to an equilibrium state 
very close to thermal equilibrium, despite the integrabil- 
ity of the problem. 

Some indications of thermalization in ID bosonic sys- 
tems have been obtained in numerical simulations of 
various physical processes in quasicondensates, such as 
the subexponential decay of coherence between coher- 
ently split quasicondensates Q, soliton formation in a 
ID bosonic system in the course of (quasi) condensation 
Q, in- trap density fluctuations @, wave chaos [l(|, and 
condensate formation after the addition of a dimple to 
a weak harmonic longitudinal confinement of a ID ul- 
tracold atomic gas However, a systematic study 
of thermalization of the GPE solution in the course of 
time evolution was lacking up to now. Even (l2| . where 
thermalization of the GPE solution with the initial con- 
ditions corresponding to the high-temperature limit has 
been numerically obtained, states that formal and sys- 
tematic understanding of the problem is still incomplete. 
We fill this gap, at least to a certain extent, with our 
present study. 

We also have to draw a clear distinction between our 
approach and that of Rigol et al. (l3j ]. who theoreti- 
cally studied dephasing in a quantum system of hard- 
core bosons on a lattice, prepared initially in a coherent 
superposition of eigenstates, and its relaxation to a gener- 
alized Gibbs (fully constrained) equilibrium. Our aim is 
to demonstrate that a weakly interacting ID degenerate 
bosonic gas can approach, in the course of its evolution, a 
state that is reasonably close to the conventional thermal 
Bose-Einstein equilibrium. 



II. NUMERICAL APPROACH 



We solve the GPE 
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where ^(x^t) is a classical complex field representing a 
quasicondensate of atoms with mass m and g is the effec- 
tive coupling constant in one dimension (we assume g > 
0). The interaction strength is characterized by the Lieb- 
Liniger parameter [l| 7 = mgj (h 2 n) = (n£)~ 2 , where £ is 
the quasicondensate healing length and n = (\^(x, t)\ 2 ) 
is the mean ID number density. We consider the weak 
interaction limit 7 < 1. We assume periodic bound- 
ary conditions for \I/(x, i), with the period L being long 
enough to ensure the loss of correlations over the half 
period: (^f*(x,t)^(x + L/2,t)) < n. The angle brackets 
denote here averaging over the ensemble of realizations. 
For each realization the initial conditions are prepared 
in a manner similar to the truncated Wigner approach 
[3] but taking into account thermal fluctuations only 
(cf. Ref. Q). We express the macroscopic order param- 
eters in terms of the phase <fi and density Sn fluctuations: 
* = (n + Sn) 1 / 2 e i<t '. The initial (at t = 0) fluctuations 
are expanded into plane waves as 

Sn(x,0) = 2 y / n/Ly^^l3 k y / rjk/e k cos(kx + zu k ), 
cf>(x,0) = (l/VnL)^2f3 k ^e k /T lk sm(kx + w k ), (2) 



where e k = yJ rj k (r] k + 2gn) is the energy of the ele- 
mentary (Bogoliubov) excitation with the momentum hk 
and r\ k = (hk) 2 / '(2m) . The real numbers /3 k and w k 
have the meaning of the scaled amplitude and the off- 
set of the thermally excited elementary wave with the 
momentum hk at t = 0. The values of w k are taken 
as (pseudo)random numbers uniformly distributed be- 
tween and 2ir. Each ensemble of realizations is also 
characterized by a distribution of the (3 k values with 
(/?^) being equal to the main number Ao(fc) of elemen- 
tary excitation quanta (quasiparticlcs) in the given mode 
|15j . In equilibrium at the temperature T the popula- 
tions of the bosonic quasiparticle modes are given by 
A/be(£,T) = {cxp[e fc /(fc B T)] - I}" 1 . 

The use of the classical field (GPE) approach is justi- 
fied, as it has been shown [l6| that the noise and corre- 
lations in an atomic quasicondensate are dominated by 
thermal (classical) fluctuations under experimentally fea- 
sible conditions, and the observation of quantum noise 
is a challenging task that can be solved in particular 
regimes by means of involved experimental tools [l7j . 

To integrate Eq. (Q}, we used the fourth-order time- 
splitting Fourier spectral method [l9j . which is rather 
similar to that used in Ref. Q- 

We have found a set of examples of solutions of Eq. 
([T]) that demonstrate quite a good degree of thcrmaliza- 
tion. Efficient thermalization has been observed in the 
cases of initial population of Bogoliubov modes within 
a certain momentum band around k = [for simplic- 
ity, we assume Ao(fc) = Ao(— k)], with the bandwidth 
being narrow enough to ensure the phononic nature of 
these excitations, |fc|£ < 1. In Fig. [TJ we present our re- 
sults of numerical integration of Eq. (TTJ) for the initial 



conditions corresponding to the truncated classical dis- 
tribution, parametrized by the effective temperature To 
and the cutoff momentum hkg, i.e., for Afo(k) being equal 
to k^To/ek for |fc| < fco and zero otherwise. For the sake 
of convenience, in Fig. Q] we plot the mean energy per 
mode E k = e k J\f(k), which does not diverge at k 0, in 
contrast to the time-dependent population distribution 
Af(k). Practically, E k can be calculated by averaging 
over the ensemble of realizations the energy stored in the 
given mode: 
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(3) 



where Sn k and v k are the Fourier transforms of the den- 
sity Sn(x,t) and velocity v(x,t) = (h/m)d(/)/dx fluctua- 
tions. 

Elementary excitations at different momenta are found 
to be uncorrelated for all propagation times, i.e., 
(5n k ,8n* k ) = (\5n k \ 2 )8 kk , and (v k >v* k ) = (K| 2 )4fc', as 
expected for a thermal equilibrium state. 

The energy distribution approaches its equilibrium, 
which is quite close to the thermal Bose-Einstein dis- 
tribution. The main difference is that the former is flat 
at k — > and the latter has a cusp there. The equiva- 
lent temperature T eq of the corresponding Bose-Einstein 
thermal distribution is determined from the energy con- 
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FIG. 1: (Color online) Dots: mean energy per mode (scaled 
to kuTo with ksTo — 2gn) as the function of wavemimber 
k (scaled to £) for the dimensionless time gnt/h, = (a) 0, 
(b) 50, (c) 2850, and (d) 5750. The Lieb-Liniger parameter 
7 = 5 x 10~ 3 , fco£ = 0.33. Solid line: mean energy per mode 
ekM[k, T cq ) for the equilibrium state, T oq = 0.35 To; see Eq. 
@. The data are averaged over 200 realizations. Units on 
the axes in this figure and the subsequent figures are dimen- 
sionless. 
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servation [181 ] : 

^e fe AA (fc) = 5> fe A/B E (fc,T cq ). (4) 

Note that for a weakly interacting ID system of 87 Rb 
atoms with the parameters as in Fig. [1] the time unit 
h/(gh) w 0.1 ms. 

To check our numerical method, we performed the fol- 
lowing tests. First, we checked the isospectrality of the 
(generalized) Lax operator of the inverse scattering prob- 
lem 0, Q . We calculated the spectrum of the linear dif- 

/ %() I Qx Q \ 

ferential operator [ ^ ^a/a^, )j where x = x/£ 



q* —id/dx 

and q = n^ 1 / 2 '^!(x 1 1), by substituting the numerically 
obtained solution for ty(x,t) at different times and com- 
paring the result to the spectrum that corresponds to 
the initial condition ^(x, 0). The spectrum of the Lax 
operator has been found to be time independent with a 
high accuracy. The maximum relative shift of an eigen- 
value over more than 100 realizations was about 10~ 7 for 
a numerical grid consisting of 1024 points in x. 

Then we checked the time independence of the numer- 
ical values of the integrals of motion of Eq. (JXJ) . The first 
three of them are (up to a numerical factor) the particle 
number, the total momentum, and the total energy of 
the system. Other integrals of motion can be calculated 
using the recurrent formula Q. We found that they are 
conserved with high accuracy, with the relative error be- 
ing of order of 10 -11 for the first integral of motion (the 
number of particles) and of order of 10~ 4 for the 15th 
integral of motion. 
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FIG. 2: (Color online) Numerically obtained energy-weighted 
squared deviation of the quasiparticle distribution from the 
Bose-Einstein thermal equilibrium as a function of time. The 
initial energy distribution and other parameters are the same 
as in Fig. QJa). The inset shows the numerically calculated 
first-order correlation function gi(x — x') (shown on the log- 
arithmic scale) for the dimensionless time gnt/h — (circles) 
and 6000 (squares). The distance is scaled to A oq . 



Following Ref. [lj|, we estimated the nu- 
merical error through the fidelity, defined as 



T = 



where 



1- {nL)- 1 J c I" dx^*(x,0)^{ h (x,t 7 -t) 

^fb(x,t, —t) is the numerical solution of the GPE with 
the initial condition ^(x,0) first propagated forward 
in time (up to time t) and then propagated backward 
over the same time interval. We obtained J- ~ 10~ 8 for 
the propagation times t as long as 10 3 Ti/{gn) 1 which is 
sufficient for the establishment of equilibrium, with the 
spatial grid consisting of 512 points. 

We found that our method converges if the grid con- 
tains more than 200 points for L w 400 £. A coarse grid 
(about 100 points) yields a numerical artifact: any initial 
distributions rapidly smears out to the "classical-like" 
flat distribution of the energy over modes, i.e. to Ek ~ 
const for all momenta — -r- < k < resolvable by the 
grid with the step Ax. 

To quantify relaxation of the system toward its equi- 
librium, we introduce the measure 



W 



E^o {e fc [A^(fc)-A/BE(fc,T eq )]} 



Efc^ [ e fc- /V BE(fc,T cq )] i 



(5) 



which has a meaning of the normalized energy-weighted 
squared deviation of the quasiparticle distribution from 
the Bose-Einstein thermal equilibrium. For the parame- 
ters of Figs. [T]and [2j with T = 150 nK, the thermal- 
ization time is r oq ~ 20 ms. If we change To to 50 nK 
and fco£ to 1, then r oq decreases by an order of magni- 
tude. Note, that the obtained thcrmalization time r oq is 
always shorter than the time needed for a sound wave to 
traverse the distance L. Therefore the thermalization ob- 
served in our simulations is a local physical effect, which 
is not related to specific boundary conditions. The ther- 
malization time T cq should not be confused with the time 
Td ~ m\\/h Q, where \t = 2fi 2 n / (mk^T) , of dephas- 
ing between two ID quasicondensatcs initially prepared 
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FIG. 3: (Color online) Dots: mean energy per mode (scaled 
to ksTo with JcbTo — 0.66 gn) as a function of wave number k 
(scaled to £) for the dimensionless time gnt/h — 2x 10 4 , which 
is long enough to provide equilibration. The Lieb-Liniger pa- 
rameter 7 = 5 x 10~ 3 , = 1.0, — 2.0. Note the close- 
ness of the equilibrium states to the initial energy distribution 
(dashed line). 
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FIG. 4: (color online) Distance D^ 2 '[t/ji,^}2] (on the logarith- 
mic scale, dimensionless) as a function of scaled time for the 
parameters of Fig. [1] (I, upper curve) and Fig. [3] (II, lower 
curve) . 



in thermal-like states with strongly mutually correlated 
fluctuations. 



III. DISCUSSION AND CONCLUSIONS 

Therefore we found numerically an example of the 
GPE solution that relaxes toward a state with practically 
measurable noise and correlation properties [20 ] well de- 
scribable by a thermal Bose-Einstein ensemble of quasi- 
particles. As an illustration, in the inset in Fig. [5] we plot 
the numerically calculated first-order correlation function 
gi(x-x') = (^*(x',t)^f(x,t))/n for t = and for t large 
enough to provide equilibration [2l|. We see that this 
correlation function finally approaches the exponential 
form gi(x — x') = exp(—\x — x'\/Xt), predicted for the 
thermal equilibrium [22j, with T ss T eq [the distance in 
the inset to Fig. [3]is scaled to A cq = 2h 2 n/(mkBT cq )]. 

Not every initial distribution relaxes toward the Bosc- 
Einstcin thermal equilibrium. For example, if there are 
initially two oppositely propagating bunches of particle- 
like elementary excitations well separated in the momen- 
tum space, an equilibrium state very far from A<be(&, ^eq) 
is established, as seen from Fig. [3J where we assume 
£fcA<o(fc) to be equal to fceTo for &i < 1^1 < ^2 and zero 



otherwise (k\ > This behavior can be viewed as a 

conspicuous example of relaxation toward the fully con- 
strained equilibrium (l3| in the weakly interacting case. 

To elucidate the qualitative difference between the 
cases shown in Figs. [1] and [3J we calculate 
the time dependence of the distance D^[-ipi, ip 2 ] = 
(2fLL) _1 Jq dx \ipi{x,t) — ip2(x, t)\ 2 between two solutions 
ipi,ip2 of the GPE, which are very close at t = 0. As we 
can see from Fig. [4] if phononic modes are initially pop- 
ulated, grows exponentially and saturates at the 
unity level (corresponding to the total loss of correla- 
tions at t — > oo), thus signifying the chaotic regime. If 
only particle-like modes are initially populated, then 
grows very slowly and stays well below 1 at all experi- 
mentally relevant times (hence, the chaotic behavior is 
practically not observed in that case). 

To conclude, we numerically observed thcrmalization 
in a ID quasicondensate, i.e. in an ultracold atomic 
system described by the NLSE with a cubic repulsive 
nonlincarity, if only phononic modes are populated ini- 
tially. The correctness of the numerical solution has been 
checked via the criteria of the Lax operator isospectrality, 
conservation of the integrals of motion, and fidelity Such 
a series of tests prevents the possible numerical artifacts 
that may occur in the split-step method [23j . Although 
the thermalization is not complete, experimentally mea- 
surable correlations are expected to be well described 
by the thermal equilibrium of bosonic elementary excita- 
tions. Our findings are in good agreement with the high 
efficiency of the evaporative cooling of ultracold atomic 
gases on the atom chips deeply in the ID regime [J, [f| 
(our work on numerical modeling of evaporative cool- 
ing of ultracold bosonic atoms in elongated traps is in 
progress). On the other hand, to provide full thermaliza- 
tion of nonequilibrium ensembles of particle-like excita- 
tions, like the one displayed in Fig. [3l we have to resort 
to the option of the intcgrability breakdown provided by 
the mechanism of effective three-body elastic collisions in 
one dimension Q. 
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